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NOTATION 


E 3 Euclidean space 

M Abstract manifold 

x k Cartesian coordinates in E 3 

v Vector in E 3 

jc" Coordinates in an open patch of M 

T p (M) Tangent space at a point p eM 

T* (M) Cotangent space at a point p eM 

V a Components of an element of T p (M ), contravariant vector 

V a Components of an element of T* (M), covariant vector 

r p ap Christoffel symbols of the second kind 

a/3 } Christoffel symbols of the first kind 
R uavfj Riemann curvature tensor 

(a i ) Hyperplane spanned by the vectors {a t } 

(a | b) Inner product of vectors a and b (abstract notation) 

u • v Dot product of two vectors in E 3 

u xv Cross product of two vectors in E 3 

f Total derivative with respect to the independent parameter, 

df / dt, df IdX , etc., depending on context 
/' Spatial derivative when/depends on only one space variable, df jdz 

V Gradient operator in E 3 

(v) Column vector with components 8 t 

(v) T Row vector, {d x d y d z ) 

D p Covariant derivative operator in M 

d 2 Total directional derivative along A , d t + A • V 

d Total derivative in E 3 , d t + r • V 

A ® B Direct product tensor of two vectors with components A j B i 

S jk Fluid shear-compression tensor d j v k + d k Vj 

co jk Fluid vorticity djV k - d k Uj 

Fatin indices, k, run from 1 to 3, and Greek indices, a , run from 0 to 3, with “0” 
reserved for the time coordinate. The Einstein summation convention is used with sum 
implied when any index appears only twice as a subscript in E 3 , and once as a subscript 
and once as a super script in M, i.e. 

A k B k = Y j A k B k , A a B a = Y j A a B a . 
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APPLICATION OF DIFFERENTIAL GEOMETRY TO ACOUSTICS: 
DEVELOPMENT OF A GENERALIZED PARAXIAL RAY-TRACE 
PROCEDURE FROM GEODESIC DEVIATION 


1. INTRODUCTION 

When background fluid motion is subsonic the acoustic rays in a time-dependent 
moving ideal fluid are equivalent to the null geodesics of a pseudo-Riemannian 
manifold 1,2 . This connection suggests a novel approach to dealing with linear acoustics. 
Equating the features of ray theory with concepts from differential geometry offers an 
advantage to traditional views by: 

• Providing a unified approach to problems in acoustics that encapsulates all known 
results in an easily accessible format. 

• Leading to immediate generalizations of known problems that include the effects 
of random media on the acoustic field. 

Specifically, by using geometric principles this approach leads to generalized 
versions of Fermat’s principle, ray integrals for layered moving media, Snell’s law, and 
the laws of geometric spreading. In addition to reproducing known results, the 
application of differential geometry leads to a generalized version of paraxial ray 
tracing 3,4 derived from the law governing geodesic deviation 5 ' 9 . The differential 
geometric approach to ray tracing also offers a new way to look at different presentations 
of ray theory as being related by a type of geometric invariance called conformal 
invariance 10 . Conformal invariance allows one to view traditional approaches to ray 
theory as belonging to a continuum of equivalent representations. 

Paraxial ray-trace procedures calculate ray paths within a ray bundle by solving a 
system of second-order linear differential equations along one ray path rather than 
repeatedly solving the ray equation with different initial conditions 3,4, '. From the results 
of this procedure the neighboring rays close to a specific solution of the ray equations 
may be mapped thus allowing one to model the deformation of a ray bundle as it 
propagates through a random environment. This type of ray tracing has become widely 
used in seismology where motion of the medium can be ignored. 

The approach presented here finds its natural home in the time domain treating 
acoustic rays as trajectories, or histories, in space and time in contrast to the traditional 
view which treats the rays as unit-speed curves in space. In the conventional approach a 
ray trace is viewed as static structure in space, sometimes called a skeleton upon which 
the field flesh lives. In contrast, the space-time approach views the rays as encoding 
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information about the history and future evolution of the acoustic field as well as its 
shape at any given instant. This information is contained in the geometry of the “acoustic 
metric” and can be thought of as a property of the medium itself for any source placement 
and geometry. 

Rays emerge as a high-frequency approximation to the wave theory of light or sound. 
The theory of bicharacteristics is more general, describing the propagation of a surface of 
initial data through space-time 22,23 . The geodesic structure that determines the rays in the 
high-frequency limit is identical in form to that of the bicharacteristics of the field 
equations. 

The rate of separation of neighboring geodesics from a specified geodesic with similar 
initial conditions is determined by the stationary values of the variation of the geodesic 
path with respect to initial conditions. This variation is measured by the Jacobi field of 
the corresponding geodesic equation determined by the second variation of the action. 
The same concept gives rise to the notion of geodesic deviation, a technique used in 
general relativity to describe gravitational tidal forces 5 ' 7 . 

Most of the applications presented in this report make use of the Jacobi field to 
directly calculate or estimate relevant aspects of an acoustic ray skeleton and the acoustic 
field. Once the Jacobi field is known for a given situation, either analytically or 
numerically, it can be used to calculate the: 

1. Field intensity in the limit of geometric acoustics. 

2. Coordinates of rays near a pre-chosen central ray within the ray skeleton. 

3. Exact location of caustics along the central ray. 

The behavior of the Jacobi field is determined by the Riemann curvature tensor 
derived from the acoustic metric—see Appendix A for details. Briefly, a positive 
curvature causes convergence of neighboring geodesics with a common point source and 
predicts the existence of conjugate points, equivalent to caustic formation in the acoustic 
field. Negative curvature leads to divergence of geodesics with a common point source, 
at roughly an exponential rate, while zero curvature leads to linear divergence. For the 
special case of constant curvature the manifolds described by the previous cases 
correspond to a sphere, pseudosphere, and Euclidean space, respectively. The relation 
between the environment and the Riemann curvature tensor allows one to determine, by 
inspection in many cases, whether or not convergence is eminent, and even estimate the 
frequency of caustic formation without explicitly solving for the Jacobi field. For 
example, stationary fluids with depth-dependent sound speeds satisfying c"(z ) > 0, such 
as the canonical Munk profile used to describe the deep underwater sound channel, are 
known to produce natural sound ducts giving rise to nontrivial caustic structure in the 
acoustic field. The curvature in such cases is proportional to c"(z), corresponding to 
locally spherical geometry. The table below provides a list of terms from ray theory and 
the corresponding, or relevant, differential geometric quantities. Each ray quantity is 
either identical to, or is derived from, the corresponding geometric quantity. 
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Ray Theory 

Differential Geometry 

Ray 

Geodesic 

Derivative of ray coordinates 
with respect to initial conditions 

Jacobi field, Geodesic deviation 

Rays and wavefronts for a point source 

Geodesic polar mapping 

Caustics 

Conjugate points 

Energy conservation, 

Fermat’s principle 

Null constraint 

Ray integrals for layered media, 

Snell’s law 

Isometry, 

Symmetry of the metric tensor 

Focusing effects of the medium 

Riemann curvature tensor 


Table 1. Ray-theory and differential-geometry analogues 

The main purpose of this report is to present a complete set of paraxial ray equations 
for three- and four-dimensional acoustic ray tracing in an arbitrary environment, 
including fluid motion and time dependence 12 ' 21 , along with illustrative examples 
employing the techniques of differential geometry, and to discuss the conformally- 
equivalent representations of ray theory. 

In Section 2 the complete paraxial ray-tracing procedure is presented in four¬ 
dimensional space-time and three-dimensional space, along with detailed descriptions of 
the quantities involved. The special case of layered media leads to an exact solution for 
the Jacobi field. Section 3 describes the effects of conformal invariance in stationary 
fluid systems and compares the results for four different popular choices of ray 
parameter. Section 4 presents an analysis of several systems that have exact solutions for 
the rays and the Jacobi field, and describes the caustic structure as predicted by the Jacobi 
field. Appendix A provides a list of relevant concepts from differential geometry as a 
quick reference for the reader. The reader interested in a deeper presentation should 
consult Refs. 7, 8, and 24. Appendix B presents a complete derivation of the geodesic 
structure associated with the bicharacteristics beginning from first principles, comparing 
the general case to that of linear acoustics in the high-frequency limit. Appendix C 
presents the technical details behind constructing a ray-centered basis for the paraxial 
ray-trace procedure. Appendix D lists the components of the Christoffel symbols of the 
second kind and the Riemann curvature tensor in a Cartesian coordinate basis, while 
Appendix E introduces isometry and derives the ray integrals along with Snell’s law for 
layered media with linear and polar symmetries. While Appendices A, B, and E are 
mostly for completeness and rigor, Appendices C and D are an integral part of the 
paraxial ray-trace procedure. 

2. GENERALIZED PARAXIAL RAY TRACING 

In this section the complete generalized paraxial ray-tracing procedure is presented in 
affine and time parameterizations along with results for layered media. Figure 1 
illustrates the basic features of a ray system in a two-dimensional slice of space-time. 
The fiducial ray is labeled y F and neighboring rays are labeled y N . The coordinates of 
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each ray are x F and x(' respectively, and are related by = x F + Y M . The choice of ray 
parameter used in the geometric approach corresponds to an affine parameter, labeled A. 
The significance of A is discussed in Appendix A.4. The ray tangent, T M , in space-time 
and the deviation vector, Y M , at two values of the affine parameter, A , are shown. The 
deviation vector lies in a 2-dimensional space-like hyperplane at each point of the ray y F . 
The Jacobi equation describes the evolution of this vector projected into this hyperplane 
and requires the introduction of a ray-centered orthonormal basis (e 7 ), I = 1, 2 (not 
shown in the figure). In general, the deviation vector will turn into the time direction 
while always remaining space-like. 



Figure 1. A fiducial geodesic and one 
of its neighbors are shown, along with 
the fiducial tangent vector and the 4-D 
deviation vector at two values of A, 
labeled 0 and 1. 


Figure 2. Similar to Fig. 1 but including the equal¬ 
time deviation vector, shown at both points along the 
fiducial ray. The 4-D and 3-D deviation vectors are 
initially equal by choice. As the system is tracked 
along y F , Y picks up a time-like component while Y 
contains only space-like components at all times. 


Coordinate time is a more natural choice of ray parameterization from the physical 
point of view. In the coordinate-time-parameterized version of the problem, the deviation 
vector remains tangent to the wavefront in space at all points along the ray—see Fig. 2. 
In this case, the deviation vector approximates a small arc in the eikonal surface and may 
be used to reconstruct the wavefront in a small region about the central ray by sweeping 
the deviation vector about the central ray for a complete set of initial launch angles, 
forming a small disk (Fig. 3). 
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Figure 3. A 3-D ray path is shown along with a small ring formed by the locus of initial 3-D deviation 
vectors. As the deviation vectors are propagated along the ray this ring deforms as a result of the local 
refractive index and fluid velocity. A pair of axes, labeled a and b, are included to demonstrate the change 
of orientation of the wavefront. 


2,1 Affine Representation 

The paraxial equations consist of three main ingredients: the geodesic equation, the 
parallel transport equation, and the Jacobi equation, respectively listed below: 


d 2 x M M dx a dx p 
dA 2 + aP dX dX 


( 1 ) 


de^ , dx a p 

-dx* x «-dx e ! =Q ' 


d% 

dX 2 


where 


r, =^g a pr p , 

K u m R m eX a X p efe V j . 


( 2 ) 

(3) 


(4) 

(5) 


The explicit forms of T^ a p and R Mav/j appear in Appendix D. Many of the longer 
equations will be found in Appendices B-D along with derivations. 
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The paraxial ray-tracking procedure may be outlined as follows: 

1. Find a specific solution to the geodesic equation, Eq. (1), labeled the fiducial 
geodesic, y F (A). 

2. Solve the parallel transport equation, Eq. (2), for the ray-centered basis e“ along 

y F W- 

3. Using the results of steps 1 and 2 solve the Jacobi equation, Eq. (3), for some 
specified initial conditions [a] to get the components of the deviation vector Y f 
along y F (A) projected onto the internal basis. 

4. The space-time components of the deviation vector, Y ft =Y I ef, give the 

coordinates of a neighboring ray, = x F + Y M . 

Expanding upon the above procedure provides the necessary information for calculating 
the geometric spread of the acoustic beam in space-time. 

5. Repeat step 3 by either: 

a. Choosing a new set of initial conditions {c?} on Y ,, such that 
Y 1 (cc,A)-Y I (a,A) = 0. 

b. Rotating the initial deviation Y(a) in the (e I ) plane by a small change of 
initial conditions. 

6. The infinitesimal cross-sectional area of the ray bundle in the (e 7 ) plane is then 
respectively given by: 

a. Y(a)xY(a). 

b. Y(a)xdY(a)/ 2. 

The order of the geodesic and deviation equations is reduced by introducing new 
variables p a = x a and Pj = Yj and defining a 20-dimensional dynamic vector, 

u = x a ®p a ®e° ®Yj ®Pj , 7=1,2. (6) 

The complete paraxial equations are then 

du EV 

—— = F(u,u) 
dA 

6 


( 7 ) 



with 


F(u,u) = p a ®-T% v p^p v ®-r Vp"c; ®Pj © -KjjYj . (8) 

The 20 equations are subject to 8 constraints. Provided the initial data satisfy these 
constraints, the solution will always satisfy them in principle. There are 7 constraints on 
the 8 components of e“ (see Appendix C) and one constraint on p a . Implementing all 
of the constraints reduces the dynamical system to 12 degrees of freedom: 

u = x k © p k © t © R{a) © T 7 © P l , (9) 

where t is coordinate time. The corresponding derivative vector is of the form 

F(u,u) = p k ®A k ®w®S>®P, © -K^, (10) 

where w = dt/dX is determined by the null constraint, Appendix B.4, Eq. (B.15), 
A k = -T k a pp a p p , with all occurrences of p° replaced using the null constraint and that 
R(a ) is an SO(2) rotation. Imposing the constraints on the system does not necessarily 
lead to a simpler system. The reduction of variables leads to an increased complexity in 
the remaining equations, and the quantities a> and w typically contain denominators that 
may vanish on occasion. 

An important consequence of this construction is that the deviation vector points in 
four-dimensional space-time from the fiducial geodesic, y F (X), to a point on a 
neighboring geodesic, y N (X ), with the same value of affine parameter X . This means 
that the deviation vector will not necessarily remain “in phase” with the fiducial ray path 
in the traditional sense of the term. This does not pose any problems in the geometric 
description of neighboring curves since the Jacobi field simply maps out the local space- 
time geometry in a tube surrounding y F (X) without prejudice to any coordinate, time 
being just another coordinate in the four-dimensional geometric paradigm. Two 
adjustments can be made which render the system of Eqs. (1)—(5) useful for 
reconstructing the wavefronts and predicting geometric loss. These two adjustments are 
independent of each other and the resulting procedure is identical to that outlined above. 

2,2 Time Parameterization 

Including the auxiliary basis described in Appendix C to the above procedure allows 
for the identification of neighboring points on the wavefront at a given value of X by 
using the data for p* 1 and e? to construct e, (see Appendix C.2, Eq. (C.3), for an 
explicit form of e r ). Switching to a coordinate-time parameterization of the problem 
allows e, to be tracked directly without the need for ej, eliminates the need for 
calculating time, and allows one to easily construct wavefronts from ray data by 
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specifying that the ray equations be integrated to a given set of time values. The 
corresponding equations and dynamical vector are listed below: 


x k = ? Sjj (x, - o i )(i ; - Vj )+ (x k - u k - In c + 2r 

•Vine) 


+ d t u k -cd k c- (o kj (i, ~Vj)+ ] 2 (s ¥ - (o kj ^, 


(11) 

—e,+ —£ iy (v x v + s x n + 2Vcx h).e r , = 0, 
dt r ’ k 2 jkA h L ' 


(12) 

£L.- K *L. + K a y l= 0, 

dt 2 dt u 1 


(13) 

with K = T°mnt m t n +r°oo + 2Y°mot” and K u ={R kl<iv 0 +R /imvn t m 

t n +2R M0m t n t 

^e v j. 

Derivations of each formula are given in the appendices. The time-parameterized 
representation of the equations requires a 16-dimensional dynamical array 

u=x©p®e I ©T 7 ©P 7 , 


(14) 

with derivative vector given by 



F = p © A ® Qe 7 ©P 7 ® -KjjYj . 


(15) 


Imposing constraints in the auxiliary basis leads to an 11-dimensional system. This is 
not a huge improvement to the 12-dimensional array of the last section. However there 
are a few attractive features of this representation of the problem. First, there is no 
equation for determining travel time along the ray path since time is the parameter used 
as a step in the integration of the system. Second, the equations are free of any divisions 
that might go to zero during the procedure. Finally, there is the conceptual or 
pedagogical advantage of being able to immediately identify the wavefronts from the data 
without any additional work. The equations do involve a greater number of terms as 
compared to the affine counterpart, introducing greater complexity. The vector T 7 
measures the deviation within the wavefront. 

2.3 Two Spatial Dimensions 

Many systems of interest in acoustic modeling are either legitimately two-dimensional 
or can be reduced to effective two-dimensional systems that provide good 
approximations to the true system. In such cases the auxiliary basis vector is completely 
determined at all points along the ray by a 90-degree rotation of the wavefront normal. 
The reduction of spatial dimensions also reduces the number of components of the 
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deviation vector from 2 to 1. In the affine parameterization there are 8 degrees of 
freedom: 

u = x a ® p a ®Y®P. (16) 

Moving to a time parameterization further reduces the system to 6 dimensions: 

u = x® p®Y ® P. (17) 

Not only is the work required to produce the answer significantly reduced, but also the 
equations are usually simpler in form for lower-dimensional systems. The basic steps in 
the procedure are identical to those outlined for the affine case. 

2.4 Initial Conditions, Mapping of Neighboring Rays, and Beam Deformation 

The components of the equal-time deviation vector Y r approximate a small element of 
arc length within the eikonal surface in the e, direction. The initial value of the 
deviation may be chosen arbitrarily. The initial value of Y r may also be chosen 
arbitrarily and is related to the values of the physical parameters that describe the 
environment. In modeling the evolution of a wavefront the initial geometry is assumed 
known, hence it is desirable to express Y I0 in terms of known quantities. By definition 
x^ = x$ + Y M , Y M = x^ - x %, and dx F / dA = dt F / dA(c F n F + v F ), with a similar 
expression for the neighboring ray, c F = c(x F ) being a shorthand notation for the local 
sound speed evaluated at a point on the fiduciary ray (with similar expressions for all 
other quantities evaluated along the ray). At the initial wavefront t F \ =t N | by choice, 
and the initial rates may chosen such that dt F /d/ 1| = dt N /dA\ . Denoting this common 
initial rate of the time coordinate by the constant (3 , the initial velocity of the deviation 
vector becomes Y 0 = fi{c N h N -c F h F +v N -v F ), all quantities being evaluated at the 
same A but along different rays. In terms of time parameterization 
dY /dt\ Q = c N n N - c F fi F +u N -0 F , and the initial deviation velocity for two rays with a 

common initial position dY / dt\ =c 0 An, where A n = n N - h F with both term s being 

evaluated at the same point in space. The quantity An determines the initial shape of the 
wavefront, whereas all other quantities are assumed to be given (see Fig. 4). 
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Figure 4. The shape of a small patch of the initial wavefront and two 
neighboring wavefront normal vectors are depicted for the three cases, 
a: (a) concave, (b) convex, and (c) planar section of the wavefront. 


The initial values of Y 1 may be expressed in terms of Y 0 M . By explicitly 
differentiating the expression Y, = Y M e lfl the following expressions for the initial values 
of Yj are derived: 

I; (18> 

where all values are taken at the initial point along y F () l) , and dY k I dA | is given in the 
previous paragraph in terms of data at the initial point, along both y F (A) and y N (A). 
Dividing by (3 in Eq. (18) gives the appropriate initial values for dY r / dt . Given a 
specific choice of initial wavefront geometry at a specific position in space-time, Eq. (18) 
gives the appropriate initial conditions for dY I / dA or d,Y I / dt. 

Intuitively caustics, or focal points, along a ray may be identified with points where 
Y M = 0. To be more explicit the deviation vector is expressed as a function of the initial 
conditions Y(A;(Y 0 ,Y 0 )) = Y A (a t ) , where (Y 0 ,Y 0 ) = a t is a shorthand notation referring to 
the set of initial conditions along the /-th neighboring ray to y F (A) (with coordinate 
indices suppressed). If the initial conditions of the family of neighboring rays are close to 
those of the fiduciary ray (which is necessarily true), then a caustic point along y F ( A ) is 
determined by Y x {a i ) = 0. It may happen that two different neighboring rays, say 
y Ni (A) and y N ^ (A) , with similar initial conditions meet without ever encountering 


d^_ 
dA 0 
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yp (A ). If the two sets of initial conditions are close to each other a crossing of this pair 
is determined by the condition Y A (a,) = Y x (a 2 ). 

By applying step 6b a picture of the beam cross section at any later time may be 
constructed. The area of the planar cross section is calculated from 

A(t) = ^Y(t)xdY(t)\, (19) 


where Y and dY may be expressed in term s of a t and da t . The acoustic intensity at 
any point along the ray is determined from the cross-sectional area of the ray bundle. 
From the differential area element dAh determined in step 5, one obtains 

IJI,=(dAn-t)J(dAn-i) 0 , (20) 

where t is a unit vector tangent to the ray path and Y is projected in the wavefront 
normal basis. 


2.5 Application to Layered Media 

To find an explicit expression for the deviation vector the tangent of the fiduciary ray 
and the basis vectors along the ray must be known at every point. For two-dimensional 
problems involving layered media it is always possible to find closed-form expressions 
for these quantities. To achieve this the explicit form of the geodesic equation, Eq. (1), is 
abandoned in favor of a set of first-order equations derived from identifying isometries of 
the metric tensor (Appendix E). 

Consider a medium with sound speed c(z) and one-dimensional fluid velocity of the 
form v x (z) = o{z). Rays, considered as curves in E 3 , that are initially fired in the x-z 
plane do not turn out of their initial osculating plane, i.e. are torsion free, and constitute 
an effective two-dimensional system. The components of the ray velocity follow 
immediately from the procedure in Appendix E: 


d't_ 
dA 


=^(i -«*) 


( 21 ) 


dx v \ 

— = K 0 a + K 0 —\i-va) 
dA c 


( 22 ) 


dz _ + rc 0 

~dA~~V 




-ua) 2 -a 2 c 2 


(23) 
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Differentiating Eq. (21) implies c 2 t = ~{k x v' + 2cc't)z, while differentiation of Eq. (23) 
and use of the last result implies z = ~{k x v' + cc'i)i , where dot denotes differentiation 
with respect to A,. It may be demonstrated that turning points in time are not allowed, as 
the relation dt I dA, = 0 leads to the unphysical result (dz I dA) 2 < 0, hence without loss of 
generality, the condition dt I dA> 0 may be placed on any ray in the system. 

The single basis vector for the eikonal, e , is constructed by simply rotating the 
normal n = dr/dt-o by ± 90 in the x-z plane, giving the following expressions for the 
components of the second-rank antisymmetric tensor A MV , e x t -e°x = ±z/c and 
e z i -e°z = + k x /c, and the useful relation e x z - e z x = ±k {) /c , which is derived from the 
identity T l e * - T k e k , = T l e, k - T k e) . Thus the explicit form of the Christoffel symbols is 
unnecessary for determining the ray tangent and its basis. 

The Cartesian components of the covariant Riemann tensor derived from the acoustic 
metric are listed below: 


R 0x0x 



(24a) 


R 0zxz =-^-r-(2c 2 u" + u(u') 2 -2cc'v') , 

Oz*z 4c 2 V V / ) 

Rq z q~ =-^y(-4c 3 c" + 3c 2 (t>') 2 +4c 2 uu"-4cc'uu'+ u 2 (u') 2 y 

The sectional curvature along the ray derived from the Riemann tensor is 

a \ a 2 ( A 2 ct , ,/ \ a 2 „ 

K = — o \Y-av) - y\°) — j-c v [l-au)+—c , 

c c c c 

which may be reduced to 


., Kq ldfld/ii (. 

K = ° , -~r\ a c "I 1 ~ av ) 

2 cdzKcdz 


The Jacobi equation then becomes 
d 2 Y 


dsl 2 


+ KY = 0. 


(24b) 

(24c) 


(25) 


(26) 


The last term in Eq. (25) states that c" governs the focusing properties of a stationary 
medium with an inhomogeneous sound speed. A ray propagating in a region with c" > 0 
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will eventually encounter conjugate points, while rays propagating in regions where 
c" < 0 will diverge from one another. When the sound speed is constant the effects of 
fluid motion can be broken up into two terms, aiu" ~(u') 2 (a/c) 2 . The second term is 
always negative, causing ray divergence. In the first term the focusing effects are 
determined by the relative sign of a and v ”. This term will cause focusing of acoustic 
rays when a and v” have the same sign. As a consequence a constant fluid velocity has 
no effect on the stability of neighboring rays, while a fluid velocity with a constant 
gradient will interfere with the focusing caused by an inhomogeneous sound speed with 
c" > 0. 

A more interesting situation occurs when both c and u depend on depth. The 
remarks of the last paragraph are still valid with the addition of an extra term coupling the 
sound-speed gradient to the fluid-velocity gradient, -u'c'K^i/c , appearing in Eq. (25). 
This term is more difficult to interpret than the others contributing the curvature. 
Consider a simple situation in which a waveguide is created by a sound-speed profile 
with c” > 0 everywhere. Above the waveguide axis c' > 0, while below the waveguide 
axis c' < 0 . If the fluid motion is to the right and characterized by a uniform gradient, 
then i/>0. For acoustic rays fired to the right a> 0, resulting in a separation of 
neighboring rays above the waveguide axis and an enhanced convergence of rays below 
the waveguide axis. When the background fluid motion is weak and slowly varying, the 
leading-order term s in Eq. (25) are av")c 2 + a 2 c”)c , indicating that the dominant 
effects are due to the concavity of the environmental parameter values. 

Between vertical turning points the deviation equation may be converted into a 
differential equation in the variable z that in general reduces to a single integral in z. 
Changing variables to dcr = cdz leads to the following general solution for the Jacobi 
field 


Y=ic { c ‘lj$ +c - 


(27) 


with constants c\ and ci set to match the initial conditions of the deviation vector. For 
cases when z = 0 identically along the ray, one cannot perform the change of variables 
from A to z and such cases must be treated separately. Equation (27) can be integrated in 
some cases to yield a closed-form solution in term s of ordinary functions. For the case of 
layered stationary media Eq. (27) becomes 


which appears in a slightly different form in Cerveny 3 . 


(28) 
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3. CONFORMALLY-EQUIVALENT RAY SYSTEMS 


Although this technique can be applied in general, we now apply it to two- 
dimensional underwater systems with range, depth and time dependence for concreteness 
and simplicity. The acoustic line element for stationary media in four-dimensional space- 
time coordinates is -c(x,t) 2 dt 2 + dx-dx = 0. In general, rays in space will not be 
confined to a two-dimensional plane but spiral through the medium as determined by the 
gradients of the sound speed. For situations in underwater acoustics one can find either 
effective two-dimensional subsystems within the ray structure, or in many situations the 
torsion of the rays is negligible for the range of interest leading to 
- c(z, r, t) 2 dt 2 + dz 2 + dr 2 = 0, with corresponding metric tensor g /a , = -c(z,r, t) 2 ® 8 tj . 

As a result of conformal invariance we are free to pick and choose from a continuum 
of theoretically equivalent representations of the same system. Why choose a particular 
one over another? Common to many laws of physics is the fact that the law itself is 
coordinate or gauge invariant yet a specific representation must be chosen to get the 
problem solved. Choosing a particular representation wisely results in a simpler problem 
from both the analytical and numerical points of view. 

Below are presented four popular choices of parameterization of the ray equations 
along with a comment of its use in the literature. Each is presented here as an affine 
parameter leading to a new geometry. In other references this is taken to be simply a 
change in parameter devoid of any geometric significance. 

3.1 Affine Parameterization 


This is the obvious choice from a geometric point of view. It appeared once in an 
article by R.W. White 19 , introducing the connection between null geodesics and acoustics 
rays. It also appears to be identical to the parameterization used by Cerveny in seismic 
ray theory 3 . The ray equations are: 


d 2 t 

1 dt 

\ dt dc ( 

dx dc dz oc't 

b 

dX 2 c dX 

[dX~8t + [ 

dX dx dX dz ) 

d 2 z dc( 

o 

II 

d 2 x dc 

(dt Y 

—- + c — 

—V c — 


dX 1 

dz^ 

dX dx 

ydx) 


The null condition may be used to eliminate dt/dX from the equations for the spatial 
variables if desired. Also, note that the second-order equation for time may be derived by 
differentiating the null constraint. When the sound speed is explicitly time dependent one 
cannot perform the single integration of dt = ds / c unless the sound speed is a separable 
function, c(t)c(r ). The Jacobi equation for this geometry is given below in Cartesian 
coordinates: 
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(32) 


\{ d z ) 

2 d 2 c | f d: 

c V d 2 c 

dz 

dx d 2 c 1 

[U A) 

dx 2 

l) dz 2 

2 dA 

dA dxdzj 


The space-time components of the ray tangent are T“ = (c“ 2 c _1 n), which may look a 

little strange at first. This is due to the nature of, and units of, the affine parameter. One 
can check by hand that (f|r) = 0 as expected. The ray velocity in E 2 is given by 
TIT 0 = cn as expected. In ray-centered coordinates the Jacobi equation becomes 


d 2 Y 

dA 2 


4^7 = 0. 

c 3 dy 


(33) 


Conceptually this is very easy to interpret. The stability of a ray is governed by the 
concavity of the sound speed in the vicinity of the ray. In practice this can be difficult to 
implement as we typically do not know a priori c(y ) (sound-speed data are usually given 
as c(x,z )). A virtue of this parameterization is the simplicity of the equations. The 
medium was not assumed to be time independent, yet no time derivatives of the sound 
speed appear in the equations for either the space coordinates of the ray or the geometric 
spreading. 

3.2 Time Parameterization 

This most closely resembles the equations found in Landau and Lifshitz 25 . When the 
sound speed depends on only depth and range, the time integral can be performed along 
the ray path indicating that the differential equation is superfluous. Time symmetry leads 
to the constraint 


d't_ 
dA 


= k 0 . 


This also indicates that travel time is a good choice of affine parameter. By simply 
rearranging terms in the line element we can find a representation of the geometry that 
allows time to be used as an affine parameter. Let dA —> dA/c 2 = dt with k 0 = 1. The 
new space-time metric structure is - dt 2 + (dx 2 + dz 2 )/c 2 = 0. The null surface is 
equivalent to a two-dimensional Riemannian surface with a line element defined by 
dt 2 = ( dx 2 + dz 2 )/c 2 . One can think of this as a travel-time space and the proper arc 
length between two points is a measure of the travel time along the curve connecting 
them. In this representation there is no equation for time. The equations for the ray 
coordinates are 


d 2 x _ 2 dx (dx dc dz <9c^| dc 
dt 2 c dt \ dt dx dt dz) dx’ 


(34) 
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(35) 


d 2 z _ 2 dzf dx 8c ^dz <9c^| ^ 8c 

dt 2 c dt \ dt dx dt 8z) 8z 

The ray tangent in space-time is T M = (l ch) with dr/dt = ch . The ray velocity has 
unit length with respect to the acoustic metric. The Jacobi equation for this geometry has 
a particularly simple form 

+ {cV 2 c - Vc ■ Vc}t = 0 . (36) 

dt 

Some attractive features of this representation are that the sectional curvature does not 
depend on the ray momentum and time is used as a parameter. Ray code based on this 
theory does not require travel-time calculations and wavefronts can be more easily 
identified from numerical data. The ray equations are not terribly complicated as 
compared with the previous case, but the spreading equation is simpler in that it does not 
depend on any of the ray-momentum variables. The second term in Eq. (36) does lead to 
theoretically interesting changes and potential numerical effects. The curvature is a 
defining property of geometry. As a simple example consider the c-linear sound-speed 
profile as a function of depth. In the affine parameterization, this is a space of zero 
curvature, i.e. flat as a plane. In the time parameterization, it is a space of constant 
negative curvature—see Section 4 for details. Additionally, sound speeds with positive 
concavity can have regions of negative curvature according to Eq. (36). In cases where 
the ray bundles of interest are required to be parabolic this does not pose any problems. 
In general, it could lead to divergent behavior of Y. 

3.3 Range-Independent Sound Speed 

This choice of parameter is used by researchers studying ray chaos in range-dependent 
systems where the ray equations are viewed as a two-dimensional non-autonomous 
Hamiltonian system 26,27 . It is important to note that in those cases a conformal 
transformation will not make r an affine parameter. This type of system has also been 
studied in the context of instability due to time-dependent fluctuations where this 
identification is possible. 

The same trick used in the last section may be employed when the sound speed is of 
the form c(z,t). In this case, beginning from the general case we have one conserved 
quantity, dx/dX = k l . This indicates that range is already an affine parameter. 
Rearranging the null-line element gives an effective two-dimensional Lorentzian 
geometry, - dx 2 = -c 2 dt 2 + dz 2 . In this representation of the ray structure, range is the 
proper arc length of the acoustic line element. Notice that the line element is negative 
indicating that the acoustics rays in (t , z) space are identical in structure to geodesic paths 
corresponding to massive particle trajectories. 
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The ray equations are 


d 2 t + 1 dc f dtY + 2 dc dt dz _ ^ 
dx 2 cdtydx) c 8z dx dx 

d 2 z +c ^ c f dt) 2 
dx 2 dzydx) 


(37) 


(38) 


The corresponding Jacobi equation is 


d 2 Y 
dx 2 H 


c dz 2 


(39) 


In this case the curvature term is independent of the ray momentum. This follows from 
the reduction of the problem to an effectively lower-dimensional system, and the fact that 
two-dimensional manifolds only have one independent component of the Riemann 
curvature tensor. The equations take their simplest form in this representation. 

3.4 Three-Dimensional Arc Length Parameterization 

This is by far the most common presentation of the ray equations found in texts on 
acoustics. It is used in the theoretical treatment of ray geometry in E 3 as a unit-speed 
curve and is the basis for a number of early ray trace codes 28 ' 30 . 

In general, the space-time line element implies the relation cdt = dS , where S is 
defined as three-dimensional arc length of the ray. This relation also immediately hints at 
an appropriate conformal transformation leading to a geometric representation in which S 
is an affine parameter. The general form of the complete set of equations is: 


dt__ 

dS ' 


1 


d 2 x _ 1 dx f 8c dz + dc dx + 1 5c^j 1 8c 

dS 2 c dS { dz dS dx dS c dt J c dx 

d 2 z _ 1 dz f dc dz dc dx 1 1 dc 

dS 2 c dS \dz dS dx dS c dt) c dz 


(40) 

(41) 

(42) 


r + KY = 0 , 


(43) 


with 
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One can plainly see how complicated things get in general. A significant feature of the 
affine parameterized equations, and the 1 © 1 -dimensional range-parameterized 
equations, is that the curvature tensor only depends on time through c and not its 
derivatives. In contrast the arc-length-parameterized representation is quite a bit more 
complex, even for simple systems. Virtues of this representation are that the ray 
momentum is (cos 6 sin 0) and the arc length is known each step of the way. 

Sample ray traces using all four approaches are given in Figs. 5 and 6, two examples 
each for the canonical Munk profile and the c-quadratic profile. A glance at the detail of 
a focal point in the c-quadratic case, Fig. 5, illustrates the precision with which each 
method located caustics. The different polygonal shapes due to different step-size 
choices are pronounced at the turning points. Both sound speeds were expressed in 
dimensionless variables and the parameters in the Munk profile were changed for 
simplicity to £ = 0.01, d = 1.3km and d = 1.5km/s . 

An embedded Runge-Kutta procedure with extrapolation and step-size 
control/prediction was used to integrate the system in all four cases. This choice was 
made to take advantage of the explicit truncation error provided by the output of the 
embedded procedure. Error control was achieved by placing a constraint of 10 6 on all 
variables in the array, including momentum. This corresponds to an error of 
roughly - 1 jus in time and - I mm in distance. The null condition was evaluated at each 
step to gauge the relative success of regulating truncation error, rather than imposing a 
tolerance on the constraints of the system. In all cases the magnitude of the error in light 
cone placement ds = 0 was bound from above by ~ 10“ 5 , and in most cases by -10 6 . 
Caustics, surface and bottom hits, as well as potential detector interactions were 
identified by a change in sign of an appropriate quantity, i.e. caustics were located by a 
change in sign of Y. The location was obtained by the method of false position. 
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Figure 5. Detail of caustic for c = 1 + 0.01^" : (a) range parameterized and (b) arc-length parameterized. 



Figure 6. Sample ray trace for Munk profile: (a) affine parameterized and (b) time parameterized. 


4. EXAMPLES 

In this section several examples having exact solutions to the ray equations and the 
Jacobi equation in terms of ordinary functions are presented and discussed. Many of the 
cases presented here may be found in the literature 13,17,21,31 and have been included here 
to shed light on the techniques outlined in Section 2. The approach presented here makes 
explicit use of the paraxial equations, Eqs. (3) and (28), with attention paid to the 
prediction of caustic curves from zeros of the Jacobi field, and interpreting the stability of 
neighboring rays by mapping divergence and convergence zones based on the sectional 
curvature. 
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In principle one knows Y{ A). The caustic is the locus of points determined by Y = 0. 
Solving this equation gives A c (a ), the value of the ray parameter at the caustic location 
in terms of the ray initial conditions. Using the ray coordinates then gives the caustic 
curve in space x c (a), z c (a), parameterized by the initial conditions of the ray. Cusps 
are determined by finding values of a for which the velocity of the caustic curve 
vanishes: 


dx c 

da 




(45) 


A few simple cases are analyzed from the differential geometric point of view. The 
simplest cases are well known examples that provide more insight into the language of 
differential geometry. Afterwards less trivial cases are presented, with solutions to the 
caustic curves determined by the geodesic deviation equation. 

4,1 Layered Media with Smooth c(z) and u(z) 


The first two examples illustrate the differences in using X and time as affine 
parameters. Both cases are described by stationary media, and the conformal metric 
introduced earlier is used. Recall from Section 3.2 that for purely two-dimensional 
systems, the sectional curvature is K = cV 2 c-Vc-Vc. For simplicity consider rays in 
the x-z plane. Along any ray an orthonormal basis may be constructed by inspection 
consisting of the unit tangent vector along the ray t = xi + zk and two basis vectors, the 
first being obtained by rotation of the ray tangent by 90 degrees, e x = —zi + xk , the other 
being e 2 =cj. In the conformal metric, coordinate time is a proper affine parameter and 
“dot” refers to differentiation with respect to time. The basis vectors (e x ,e 2 ) are 
orthonormal with respect to the conformal inner product g tJ =8 ij lc 2 . The sectional 
curvature in the (t,e x ) plane is K = cc" - (c'f , while the sectional curvature in the 
(I, e 2 ) plane vanishes. From these results it follows that rays fired from a point source at 
the same initial angle but in different planes will diverge linearly in time. 

Spherical Geometry and Ideal Focusing 

By direct calculation the sound-speed profile c(z) = C cosh (^/x^z/c) produces a 
space of constant positive curvature K = K 0 . The deviation vector can be solved exactly 
for all rays in this system: Y(t) = A cos(VaT) + Bsin(-JlCt ). One can conclude from this 
that ideal periodic focusing will occur along each and every ray with the same frequency 
regardless of the placement of the source. The time-versus-depth integral can be 
calculated in this case leading to the result that the travel time between turning points is 
independent of initial conditions, hence the point-like focusing demonstrated in this case 
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will always occur 31 . In the affine representation the sectional curvature is K = a 2 K 0 . 
The appearance of the ray parameter leads to a nontrivial change in the interpretation of 
the geometry of space and its effect on geodesics. 

With time parameterization, the geometry is identical to a two-dimensional sphere. 
The periodic ray structure can be mapped into a sphere by morphing the rays in such a 
way that they resemble those of a source placed on the waveguide axis, then identifying 
consecutive focal points with the north and south poles of the sphere. The rest of the 
space is then projected stereographically onto the sphere. The resulting space, while 
being geometrically identical to a sphere, is topologically equivalent to a sphere with two 
holes pinched in it identified with z = ±oo . 

With affine parameterization, the geometry is spherical along each ray with a different 
effective radius of curvature for each value of a . Rays fired towards z = ±oo never turn 
and their neighbors drift away from them at a constant rate. In the previous case all rays 
experience the same curvature, indicating that even those fired towards z = ±oo converge 
with their neighbors passing through caustics. It can be easily verified by direct 
calculation of the ray integrals that the total travel time from - oo to + oo is equal to the 
period between consecutive conjugate points. Hence for a source placed anywhere in the 
space, the vertical rays meet the end of their journeys before ever encountering their 
neighbors. The finite travel time from -oo to +oo is an indication that the space 
described in the first case suffers from geodesic incompleteness 13 . 

The Pseudosphere 

The exact ray paths for the sound speed c(z) = a+bz are circular arcs. The induced 
metric space is known in the literature as a Lobaschevsky space and is an example of a 
space of constant negative curvature, or pseudosphere, K = -b 2 . The deviation vector in 
this case may be found by inspection: Y(t) = A cosh(ht) + Bs\n\\(bt) , implying that the 
distance between neighboring rays grows exponentially in time. Once again the affine 
version of the theory predicts a slightly different result. The curvature in this case is 
identically zero indicating that rays diverge linearly in A. 

Harmonic Waveguide 

The harmonic waveguide, described by n 2 =a + bz 2 , has been the subject of many 
studies due to the existence of exact solutions with nontrivial behavior for the rays and 
modes associated with this profile 30 . Consider the special case c 1 = -Jl-(l-z) 2 , with 
the units of speed scaled to 1. This sound speed forms a waveguide with symmetry axis 
at z = 1 and c —> oo at z = 0 and z = 2. The exact ray paths are sinusoidal: 


z = 1 + *J\-p 2 sin(A + <p 0 ), 


x = pA , 


z 0 — 1 



(46) 
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Consider a point source placed on the waveguide axis, z = 1. The sectional curvature is 
positive definite, the corresponding Jacobi field is 


Y = P 2 sinA + (l-/7 2 )AcosA ^ 

tJ\ - (1 -/> 2 )sin 2 A 

For the special cases p = l and p = 0, 7 = sin A and 7 = A, respectively, as expected. 
The caustic is determined by solving the transcendental equation 

tan A = A. (48) 

P 

A plot of the rays and the caustics for a source on the symmetry axis is given in Fig. 7. 
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Figure 7. Ray fan and caustic structure for the n 2 quadratic sound-speed profile. The rays are plotted for 
even increments of cos 6 0 • The caustics are traced by inserting solutions to Eq. (48) into Eq. (46). 


The next few cases involve fluid motion. Very few interesting cases with exact solutions 
exist. We consider first the curvature induced by a class of symmetric velocity profiles 
that are found to create ducted regions. These ducts are identified by mapping the 
regions of positive sectional curvature into the x-z plane. Afterwards, the exact ray paths 
are presented for two special cases with nontrivial sound speed and fluid velocity. 
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Wind-Induced Sound Ducts 


In general, effects of u' on a ray system are illustrated for a class of velocity 
functions of the form u(z) = c 0 (z/L ) 2 ”, where c 0 is the local sound speed, assumed to be 
a constant L , a length scale chosen such that the flow becomes supersonic for |z| > L and 
n > 0 is a positive integer. This class of functions has the following common properties: 
o{z) is a monotonically increasing function forz > 0, v(z) > 0 for allz, v'(z) >0 forz > 
0, o”{z) > 0 for all z, and u(-z) = v(z ). For the time being attention will be paid to rays 
that remain in regions of subsonic flow. 

From Snell’s Law, e{z) + sec# = s 0 + sec# 0 , where s 0 and # 0 are the Mach number 
and initial angle of inclination at the source point, respectively. For purposes of 
illustration it is assumed that the source is place at z = 0 and attention focused on rays 
fired in the upper half plane. For rays fired against the wind, # 0 > jt)2 and sec # 0 < -1. 
Combined with Snell’s law this information implies that sec# < 0 at all points along the 
ray, indicating that rays with a negative ray parameter will never turn back toward the 
source. On the other hand rays with positive ray parameter will always encounter vertical 
turning points, and because of the symmetry of o{z) this will happen periodically. 

The vertical turning points may be found from Snell’s law by setting 6 = 0 and 
solving for z. Although the form of Snell’s law used here is expressed in terms of the 
direction of the wavefront normal, the geometry of the rays indicates that 6 = 0 when the 
ray encounters a vertical turning point. The equation for the turning points becomes 
s(z) = -1 + sec # 0 . The condition 2 > sec Q 0 => # 0 < ;r/3 is imposed to single out rays 
that turn before entering supersonic regions, a condition that holds for any choice of v. 
The sectional curvature for this class of problems is 

K = ^^2n(4n-\)C 2n ~ 2 (49) 

L c 0 y cos 6 An-\) 

for all rays in the system, where the dimensionless variable £ = z!L is used. Clearly the 
sign of K depends on two factors, cos # 0 and the term in parenthesis. Equation (49) will 
change sign when cos 0 = cos # 0 (2 -1 / 2n ), the same result holding for rays fired both 
into and against the wind. Since cos# is an increasing function along rays with 
cos # 0 < 0, the condition is never satisfied and these rays will remain in a divergence 
zone forever. Rays with cos # 0 > 0 are a little more interesting. The curvature will 
remain positive along these rays only if cos# < cos# 0 (2-1/2 n). Since u is strictly 
increasing cos # is strictly increasing as well and it suffices to check the curvature at the 
turning point where cos # = 1. The curvature will be zero at the turning point if 
cos# 0 = 2n/(4n -\). For n = 1 this leads to the constraint cos# 0 =2/3, or # 0 »48°. 
Rays fired into the wind with cos # 0 >2/3 will remain in a region of positive curvature 
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as they propagate downstream. Rays with cos 0 O <2/3 will start with K > 0 and 
eventually along the ray K < 0, causing neighboring rays to begin diverging with respect 
to A, in general when cos 0 = 2n cos 0 () /(An- 1). Figure 8 shows the convergence and 
divergence zones for a point source placed on the waveguide axis. The zones are labeled 
by the sign of the sectional curvature, K. 

In the case n = 1, the ray equations and the deviation equation may all be integrated 
between turning points to yield solutions for the time of flight, range, and deviation 
vector as functions of depth in term s of elliptic functions. A ray trace using the same 
technique employed in Section 3 is shown below along for one half of a ray fan, with the 
placement of the caustics along each ray shown in the detail (Fig. 9). The procedure was 
stable with an upper bound on the deviations from the null surface ~10 4 for the same 
truncation-error tolerance. Figure 8 shows the regions of positive and negative curvature 
mapped in the x-z plane for a source at the origin. There is a cut-off for rays fired with 
the wind at = 2/3, z = 1 / a/2 . Rays with p > 2/3 have turning points at z < 1/ a/2 and 
remain completely within a region in positive curvature. 



Figure 8. Convergence and divergence zones for the quadratic fluid-velocity profile. Regions are 
marked according to the sign of K. The outer K = 0 boundary corresponds to the CC = 0. The inner 
boundaries correspond to a source placed at the origin, on the symmetry axis of the guide. 
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Figure 9. Ray trace for (left) the quadratic wind profile, with detail (right) for one cycle of the ray- 
path motion illustrating caustic formation. 


Linear Fluid Velocity 

For a linear fluid-velocity profile, v = czlL , c = constant, and K = -k\/ L 2 . The 
deviation equation may be integrated to yield 

Y = (t 0 /4k )sinh(VZ;t) + Y 0 cosh (JkX) . (50) 

Hence neighboring rays with the same initial position will diverge at a rate of 
Y = (f 0 /)sinh(a at, A / c) with respect to the affine parameter. The ray coordinates 
may also be solved in this case. Defining new dimensionless variables along the ray 
£ = a (eg) -1, g= z/L , t = c 2 a t/L, rj = x/L, a = aXjL , and the new ray parameter 
/? = ca , Eqs. (21)-(23) may be integrated to yield 

^ = P cosh(cr - (p) (51) 

t = /?sinh(cr-^)-^ 0 (52) 


^ = 4+ ^ ln A+^°'-^( 2 -Acosh(cr-^))sinh(f7- <p) . 


(53) 


where the parameters A 0 = ^ 2 - (3 2 +£ 0 and (p = ln(/?/T 0 ) have been defined. Eqs. 
(51)-(53) represent the affine-parameterized coordinates for a segment of the ray with 
z > 0 and the initial conditions set so that r(0) = 0. The ray-path geometry implied by 
these equations is well known 13 : 
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( 54 ) 


, = +1Lnf 

2 ^(l-0-V(l-0 2 -a 2 J 
- ^ t ((1 - a £) 7 (i - ^) 2 - « 2 - (! - a ^ o )^ 1 - a€of-a 2 ) . 

Equation (54) gives the upward part of the ray while the downward portion is determined 
by symmetry to be 2rj tp - 77 ; tj 0 , are the initial coordinates of the ray and ij tp , % tp are 
the coordinates of the turning point. 

Linear Fluid-Velocity and Sound-Speed Profiles 

Generalizing the last case, we consider the case v = £ and c = 1 + sE, . The exact ray 
paths are 


a 2 (£- 2 -1) 


7(1 -aof 


f (1 - ao) + 
H 1 . 


£ 2 sls 2 -1 y(l-au) 2 -a 2 c 2 sls 2 -1 


~U U 1 r -r-r-clIULClIl , ,— 

a (^ -1) ^ VO - a °f ~ a2cl J £l ~ 1 

(1 -au) + ca£ J 

“J 

(1 -au) + 7(1 -ao) 2 -a 2 c‘ 


H—j- In 2(# + £) — 


7 = — t - - J ( l - ao ) 2 - a 2 c 2 

L a£ k 


while the Jacobi field along the ray as a function of depth is 


^ = 0,(1 a °) +c ^_ ao f _ a i c i 

a (a + £) 


(55) 


(56) 


The sectional curvature is K = -a(a + £) /c 3 . In this case the curvature is negative for 
rays fired with the wind, a > 0, indicating divergent behavior. For rays fired against the 
wind the sign on the curvature depends on the sign of (a + £ ). Setting (a + £) = 0 
implies cos 0 O = -£ . For £ > 1 the fluid motion is always subsonic and a change in sign 
of K for rays fired against the wind is not possible, while for 1 > £ > 0 fluid flow 
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becomes supersonic at z = 1/(1 - e) , K > 0 when |cos# 0 | < e and K < 0 when |cos <9 0 | > s . 
The curve defined by a = 0 is 




M I 



and the curve defined by a = -s when 1 > s > 0 is 


la=-e 


V l-E 2 


( 1 - 4 )- 



A plot of the convergence zones for each case is given in Fig. 10. 


(57) 


(58) 



Figure 10. Convergence and divergence zones for the c-linear, V -linear case outlined above. 

(a) corresponds to £ = 0.5 and is sound-speed dominated with subsonic flow, while in case 

(b) £ = 1, with supersonic flow occurring at E, = 1. 


4,2 Piecewise-Linear Profiles 

Although the linear sound-speed and fluid-velocity profiles lead to ray systems with 
somewhat benign behavior, they are frequently used to formulate exact ray theoretical 
models of real situations where the sound-speed and velocity profiles are approximated 
by piecewise-linear functions. In such cases the ray paths and their first derivatives are 
continuous and exact solutions to the ray path and the Jacobi field are known in each 
piecewise-linear region. The discontinuity in the sound-speed and velocity gradients 
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leads to the presence of a Dirac delta function in the sectional curvature, in turn leading 
to special boundary conditions for matching the Jacobi fields of different segments at the 
boundary between two regions. The basic situation is illustrated in Fig. 11. The medium 
is divided in depth into a finite number of regions. As a ray travels through the medium 
it will turn and periodically encounter each region, in principle an indefinite number of 
times. The ray segments are indexed in order of increasing ray parameter without 
reference to the corresponding region. The Jacobi field is a function of the affine 
parameter evaluated along the ray path. Each ray segment in space corresponds to an 
increment of the affine parameter. Passing from one region to another in space 
corresponds to boundary point along the A axis—see Fig. 12. 



Figure 11. A physical situation approximated by piecewise-linear sound-speed and 
fluid-velocity profiles. A sample ray is drawn illustrating the meaning of the index 
labels along the ray as well as those labeling the different regions. 


Once the rays are matched up, the parameter value, time of flight, and range for each 
segment may be determined. The Jacobi field is continuous leading to the boundary 
condition Y i+l = Y i . Integrating Eq. (26) across a boundary leads to the following 
condition on the discontinuity in Y : 


28 





y <+1 (A a )-y < (A Jl )+0y < (A a ) = o, 


(59) 



( aAu'(l - au B ) + a 2 c B Ac') , 


(60) 


where the subscript B means evaluated at the boundary and At/, Ac' are the 
discontinuities in the fluid velocity and sound speed, respectively, at the boundary. Once 
these conditions are matched for a significant number of segments the caustics and their 
features are determined in the same manner as before. A striking feature of these models 
is that nontrivial focusing emerges as a result of the boundary conditions, exclusively the 
curvature being otherwise zero or negative. In general the rays and the Jacobi field in 
each region are given by Eqs. (55) and (56) of the last section with appropriate changes in 
scale of the coordinates. 



Figure 12. The Jacobi field for the profile illustrated in Fig. 11. The individual segments 
of Y shown here are linear merely for simplicity. The different depth regions in Fig. 11 
correspond to intervals on the A axis. Comparing Figs. 12 and 11 indicates, for 
example, that the 4 th and 6 th intervals in Fig. 12 both correspond to region 3 of Fig. 11. 


Two special cases of interest are: (1) a stationary medium with a piecewise-linear 
sound speed and (2) constant sound speed with a piecewise-linear fluid velocity. In case 
1 the curvature is identically zero in each region and the Jacobi field may be written 
Y i =A i + A, while in case 2 the curvature is a negative constant in each region and the 
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Jacobi field is Y f = A i e'^ A + B t e ' / ~ AA . In each case applying boundary conditions on Y 
leads to expressions for the coefficients A and B in terms of the ray parameter p and the 
source coordinates. The caustic in each region is determined by the condition Y i = 0 
leading to A iC , with A ;C =-A i /B i in case 1 and 2^-KA iC = ln(- B i /A i ) in case 2, for 
the value of affine parameter along the z'-th segment of the ray at which the caustic 
occurs. Inserting this expression into the appropriate ray equation then gives the caustic 
coordinates along the ray as a function of initial conditions. The piecewise-linear profiles 
give rise to cusp and broken caustics depending on the relative values of the profile 
slopes in consecutive regions. Two cases utilizing this procedure are illustrated next. 
Both cases involve stationary media due to their simplicity and past use in underwater 
acoustic modeling. The generalization to non-stationary media based on the steps 
outlined above is immediate. 

Sound Incident on a Linear Half Space 

This classic problem of sound in a homogeneous medium incident on a region of 
space with a graded refractive index is known to produce a caustic in the emerging 
field 31 , and in the graded medium. In affine parameterization, the equations for the ray 
path are 


dx dt 1 dz yl\-a 2 c 2 

— = a, —, — = ±-, (61) 

dA dA c 2 dA c 

with the local sound speed given by c(z) = c o 0(z)+ c 0 (l-z/Z)0(-z), where 0(z) is 
the Heaviside-step function and the dimensions of A are L 2 /T. Defining cos 0 O = p the 
rays in region z > 0 are straight lines given by 


rj = pA + rj o, £ = £ 0 ±jl-p 2 A, t = A + t 0 . (62) 

Here we define the dimensionless variables A = A / Lc 0 , £ = z / L, tj = x/L, and 
t = tc 0 / L . For simplicity we set r 0 = 0 and consider a source placed at r/ 0 = 0 and 
> 0. Furthermore, we focus attention on rays with sin 0 O < 0. The solution of the ray 
equations in the linear region may be expressed as 

(? - (/, + P' Vl -P ! | + (f -1) 2 = p- 1 , (63) 

where ij l = cot 0 {) . The equation for ?; is the same in both regions. The ray paths in 
space are circles of radius r = |sec <9 0 1 centered at ( r/ x + tan 0 (] l). By symmetry the ray 

path emerges from region 2 at - 0 {) . Hence the third segment satisfies Az/Ax = - tan 0 {) . 
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At the second boundary point we have rj 2 = cot 0 0 +2 tan O 0 , leading to the following 
equation for £: 


{(A) = J^]SA-({ 0 +2(p - 2 -l)) . (64) 

The initial conditions on Y 1 are chosen to describe the geometry of neighboring rays 
from a point source, Y l =80 A. Applying the boundary conditions leads to the following 
recursion relation 


Y t (A) = Y t _ x (A) - (A M )(A - A w ), £~2, 3 , (65) 

for the solution along each segment of the ray path, with k = p 2 /^/l- p 2 . Solving for 
the coefficients gives A 2 = 1 - kA 1 , A 3 = -1, B 2 = xA 2 ,, and 5 3 = 4/r 1 , with 

Y 2 = (\-kA,)A + kA], Y 3 =- A + 4k~ 1 , (66) 

for the Jacobi field in each region. An overall factor of SO has been divided out for 
convenience. For a caustic to exist in the lower portion of the medium it is sufficient to 
require K 2 (A 2 ) < 0 . This places the following constraint on the initial conditions of the 
ray path: > 2 [p 2 -l). This condition implies two things depending on the question 

asked. First, for a given initial ray labeled p, at what location of the source does the 
caustic on that ray pop out (in) the lower half space? For > 2 (p~ 2 -l), the caustic 
along ray p will be inside the lower region. Asked another way, for a fixed source 
location which rays pass through the caustic in region 1 and which in region 2? Solving 
for p in the inequality gives 



Rays fired at steep angles (small p) will penetrate deeply and pass through the caustic 
in the upper half space while those with shallow launch angles (large p) will pass through 
the caustic in the lower half space. These statements assume that the ray in question will 
pass through a caustic. 

Along ray segment 3, we have A c = 4 k 1 . Using the ray coordinates gives 

4 = 2 ( p - 1 r ,, = A />‘ 2 -1 ■ ( 68 ) 


or simply 
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( 69 ) 



for the caustic curve in the upper half space. The caustic curve in the lower portion of 
space is more difficult to extract. For a caustic to form along segment 2 of the ray path 
we have A c = kAj/(kAj -1). To extract a curve /(£,, rj c ) = 0 for the caustic, 
independent of the initial conditions, requires solving the following equation for p: 

(a 2 +b 2 )p 6 -b(2b + \)p 4 + (2b+ l)p 2 -1 = 0 , (70) 


where b = 1 + and a = £ 0 2 / rj c . The exact roots may be found by substituting p 2 =q 
and solving the third-order equation. Rather than solve for ^ c (p c ) in region 2, we 
express the caustic as a parameterized curve (^ c (p),T) c (p)) for fixed £ 0 : 


The velocity of the caustic is 

dq c _ g2 P 2 (p 2 ^ o+3)-3) 
d P (l-y) 3 / V(£ 0 +l)-l ) 2 ’ 

d^ = ^ 2 _ p{p\Z o+3)-3) _ 

dp ° (tf o+ l) p 2 -lf ,2 ( p * + p 2 (tf +2{ 0 -2) + l-2{ 0 ) 


(72) 


(73) 


Clearly both vanish when p = ^3 /(£ 0 + 3) indicating the development of a cusp. This 
cusp will always exist in region 2 as long as > 0 . Ray traces and caustics are shown 
for a few source placements in Fig. 13. The caustics, appearing on top of the ray trace, 
were derived from the solutions presented here. One can clearly see the development of 
the cusp in the lower half space. As the source is pulled up, the caustic sinks further 
down into the lower region with the cusp moving further away from the origin. The tail 
of the caustic approaches zero as rj c —> 00 . In the limiting case where the source in 
placed on the 7 axis, the cusp moves right up to the source and the tail merges with the 
7 axis. 
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Figure 13. Ray fan incident on a c-linear half space for 3 source placements, from upper left going 
clockwise, £ = 1,5, and 10. 


The Wedge Duct 

Consider c(£) = (1 + £)0(£) + (1 - o£)©(-£) (Fig. 14). The upper portion of the duct 
is fixed at a 45-degree slope, while cr measures the steepness of the lower part of the 
duct relative to that of the upper part. For simplicity consider a source placed on the 
positive £ axis. The initial ray parameter is p = cos# 0 /(l + £ 0 ). The ray is divided into 
segments labeled y i . The points where the ray crosses the ij axis are labeled rj i and the 
corresponding parameter value A,.. The first segment is given by 
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{rj + p^l-p 2 (l + £ 0 ) 2 ) +(l + ^f=p 2 . 


(74) 


with a - sign for rays fired upward and a + sign for rays fired downward. After 
traversing this segment, the ray periodically crosses the rj axis with the same value of ray 
parameter determined by Snell’s law, cos0, = cos# 0 /(l + £ 0 ). In region I the ray is 


(?7 - (?7, + ^-P 2 )J + H = p 2 <t 2 , 

while rays in region II are given by 

(7 - ( 7 , + P 1 + (I + £) 2 = P 2 ■ 


(75) 

(76) 


The crossing point of the first segment is r\ x = p l y]\ - p 2 ± p 1 ^j\ - p 2 (\ + £ 0 ) 2 , and the 
crossing points for the other segments are 


- r l\ + y (l + cr l )-2)p \ 

h~p 2 , 

7 = 2, 4, 6, 

= 77l + (7 - l)(l + 0 - 1 )p 1 V 


7 = 3, 5, 7, 


for the lower and upper arcs, respectively. Applying boundary conditions leads to the 
following expressions for the coefficients A t and B i : 

A _ J 1 ±( 1 '—1)(1 + <7 _1 )R 7 = 1,3,5,... 

' _ {- o-(l ± (7 - 1)(1 + cr- 1 )R ) 7 = 2,4, 6 , ... 

_(1 + ct )(._ 1 X 1 ±(._ 1)(1 + ct . , );; + ;( 2 ) 

s,= m , 

( +<T) ((i -1) ± (7(7 -2)+ 2 + 7(7 -2 )07? + (7 -1)7? 2 ) « - 2,4, 6 ,. .. 

I tc 


where, for simplicity, we have defined the ratio 


Vl-Ql + O 


and k = p 2 /-y/l — /? 2 as before. 
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Sample rays for cr = 5 and £ 0 =1 are shown in Fig. 15 along with caustics. The 
caustics are plotted for q e [0.2,0.9]. For q = 1 the upper and lower caustics terminate at 
ranges of 4.15 and 2.08, respectively. Both become unbound as q approaches 0. The 
Jacobi field along the first three segments is displayed in Fig. 15 for the four rays in Fig. 
16. 



Figure 14. Sound-speed profile for the wedge duct. The lower portion is plotted for cr = 5 . 



Figure 15. Caustic for the bilinear duct. Only 4 sample rays were traced for cleanliness. 
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Figure 16. The Jacobi field plotted for the first 3 segments of 4 different rays, cr = 5 and <^ 0 = 1 . 


5. DISCUSSION AND CONCLUSION 

The space-time dynamic ray-tracing procedure outlined in Section 1 provides a 
generalization to the approach currently used in seismology, and offers a method of beam 
tracing in a random medium. The identification of acoustic rays with null geodesics of a 
pseudo-Riemannian manifold leads immediately to this generalization. The geodesic 
deviation vector, used to model the deformation of an acoustic beam, allows one to 
measure geometric transmission loss and caustic location without repeatedly solving the 
ray equation. Additionally this identification has led to a continuum of theoretically 
equivalent paraxial ray-trace systems connected by conformal transformations, and points 
to a unified approach to extracting well-known results such as Snell’s law and Fermat’s 
principle in random fluid media. The null hypersurface is viewed as a manifold created 
by the local sound speed and fluid velocity which encodes all information about the 
acoustic field. 
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The development here was intended for describing acoustics in the presence of 
subsonic flow. This requirement is necessary for the identification of acoustic rays with 
null geodesics. For cases involving supersonic flow, the bicharacteristic condition 
defines the geometry of the acoustic rays as Euclidean line elements of zero length. Such 
spaces are a special case of a wider class called Finsler spaces. This identification 
indicates that these techniques may be generalized to include cases involving supersonic 
flow. 
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Appendix A 

CONCEPTS FROM DIFFERENTIAL GEOMETRY 


A.l Metric 

The metric is a generalization of the dot product and is usually introduced via a 
quadratic form 

ds 2 = ciydXjdXj , (A.l) 

with ds a differential arc length and a ij {x k ) continuously differentiable functions of 
coordinates. One can always find a small neighborhood about an arbitrary point such that 

a tJ ~ diagQ^A -1... -1) , (A.2) 

modulo constant scale factors, with n + m = dim(AZ). 

When m = 0, the metric reduces to the identity matrix and the geometry is locally 
Euclidean. When \n - m\ = dim(M) - 2, the metric reduces to one of two choices for the 
Minkowski metric used in special relativity, specifically for dim(M) = 4. Manifolds with 
metric defined by Eq. (A.l) that are locally Euclidean are referred to as Riemannian, 
while those that have local Minkowski geometry are termed either pseudo-Riemannian or 
Lorentzian. 

A consequence of Eq. (A.2) is that curves on a pseudo-Riemannian manifold can have 
positive, negative, or zero length. The coordinates of the manifold are divided into x t , i 
= l...dim(M) - 1, and t, and choosing m= 1. With this convention, curves with ds 2 > 0 
are called space-like, and curves with ds 2 < 0 time-like. These two types of curves are 
separated by a hypersurface defined by the class of curves with ds 2 = 0 referred to as a 
null hypersurface, a generalization of the light cone. 

For subsonic flow, the metric introduced in Appendix B describes a Lorentzian 
manifold with n = 3 and m = 1. The characteristic condition identifies the acoustic rays 
with geodesic curves on the null hypersurface similar to photon trajectories in general 
relativity. 
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A.2 Manifolds 


A manifold is an abstraction from the theory of surfaces that does not include a 
surrounding Euclidean space of higher dimension in its description. For a detailed 
definition of a manifold, see the References. We present here only a few concepts that 
are used in other places within the text. A tangent space, T p (M), is defined at each point 
on the manifold which serves as an abstraction of the tangent plane to a surface. Each 
tangent space has the structure of flat Euclidean space, modulo constant scale factors. 
The metric in each tangent space, g (x °), is fixed but is allowed to vary from point to 
point in M. If U,V eT p (M) then their inner product is given by g flv U fl V v , and 
g flv V fl V v = ||f|| 2 defines the magnitude of a vector. The V M are the components of the 
vector in some chosen basis of T p (M). A dual tangent space, T* (M), is also defined at 
each point of M. The metric tensor defines a mapping from T p (M) —» T*(M ), 
T =g MV T v . The T are the components of the dual vector in some chosen basis of 
Tp(M ). Tensors of arbitrary rank are elements of a direct product space constructed of 
multiple copies of the tangent and dual tangent spaces 
T p ( M ) x — T ( M ) x T p ( M ) x---xT* ( M ). For example the metric presented in Section 
B.l are components of a tensor in T* x T *. These spaces do not require a concept of 
distance to be defined. In some sense defining a metric before this section is like putting 
the cart before the horse. The components of a tensor of arbitrary rank are denoted 
T a '"' aN p l -p M . A standard choice of basis for T p (M ) and T*(M ) are d p and dx u 
respectively, where x M are coordinates defined on an open set of M containing p. In this 
basis a vector and its dual are denoted V = V^d p and V cb c ,u , respectively. 

A.3 Covariant Differentiation and Parallel Transport 

The covariant derivative of a vector and dual vector are defined by 

D fl U v ~ d p U r + , D m U v ~ d p U v - T \ V U X , (A.3) 


respectively, where 

r V ^g VP {dpgxp + dxgp P -d lj gp X ). (A.4) 

These Christoffel symbols take into account the turning of the basis vectors as one 
differentiates the vector, requiring comparison of V at two points of M. Eq. (A.3) is a 
component form of the result. 
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Comparing vectors at different points in M is not a trivial operation since they live in 
different spaces. To compare V e T p (M ) to V e T q (M ) , p,q eM with p * q , one of 
the vectors must be moved into the space of the other. In general this requires specifying 
the path in space, C, traversed by V from p to q. If T^ 1 is the velocity of C then V is 
parallel transported along C if 

T M D /J V v = 0 . (A. 5) 

A.4 Affine Parameter 

Geodesic curves parallel transport their velocity vector, i.e. 

DT V 

-= T M D T = 0 . (A.6) 

dA 

When Eq. (A.6) holds the parameter A defined along the curve is called an affine 
parameter. Any other choice of parameterization will change Eq. (A.6) to 

T f, DJ v = oT v , (A.7) 

where <7 is an arbitrary scalar function defined along the curve. Equation (A.6) makes 
the transported tangent vector “parallel” to the tangent of the curve at the new point 
(allowing a change in magnitude), whereas the auto-parallel condition requires the 
transported tangent vector to be identical to the tangent vector at the new point along the 
curve. For space-like geodesics, arc length is an affine parameter. Time-like geodesics 
may also be parameterized by their length interpreted as a standard internal time, or the 
proper time of an observer whose world line coincides with the particular time-like curve. 
Clearly, for null curves arc length cannot be used as a parameter. In such cases one 
imposes Eq. (A.6) and the null constraint thus defining the parameter as affine with no 
physical significance attached. 

A.5 Geodesics 

Defined in the previous section as auto-parallel, the geodesic equation is also the 
Euler-Lagrange equation derived from a variation of the arc length, 

S\ds = 0 , (A. 8) 

with ds defined in Appendix A.l. Equation (A. 8) defines geodesics as curves of 
“optimal” length, either minimizing or maximizing the length between two points of M. 
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A.6 Curvature 


The Riemann curvature tensor, and the covariant Riemann tensor are derived from the 
Christoffel symbols: 


R“ avf > = d v Y p a p -d p Y M av +T t ‘xvT\ p - T W, (A.9) 

R„ av p = d v {p\ap)-d p {p\av}+Y p av {p\pp)-Y p aP {p\pv). (A.10) 


The tensor defined in Eq. (A. 10) has several symmetries that reduce the total number of 
components: 

R-nvaP = Rapnv = ~^vjuap = ~^uvfla • 


The effects of curvature are most easily seen by parallel transporting a vector around a 
closed loop in M. If this task were attempted on a plane or in Euclidean space the result 
would be trivial, the transported vector would be identical to the original vector. A 
simple nontrivial example of parallel transport on a curved manifold is constructed by 
moving a vector from the north pole of a sphere around an equilateral geodesic triangle 
with 9 = 90°. The transported vector is rotated by 90 relative to the old vector. 


A.7 The Jacobi Equation 

The geodesic equation gives curves of optimal length on M satisfying Eq. (A.5). The 
stability of a geodesic relative to slight variation in initial conditions is governed by the 
Jacobi equation. For convenience we define the Jacobi vector, Y, along the geodesic 
subject to the constraint g fjv Y p x v = 0. Focusing exclusively on the null geodesics of M 
confines Y to the space-like plane (e ,). Geometrically Y points from one geodesic to a 
point on a neighboring geodesic with the same value of A . The vector Y obeys a second- 
order linear differential equation along the geodesic: 


d% 

dA 2 


+ K IJ Y J = 0, 


(A.l 1) 


which is derived from the second variation of the length integral. T 7 = e“ g ap Y p is the 
projection of Y onto the geodesic-centered frame field defined in Appendix C, and the 
curvature matrix K u = R^ avp x a x p e p e V j is introduced. For I = J, K u , no sum on /, is the 
sectional curvature in the plane spanned by (i,e 7 ), and measures the relative rate of 
separation of neighboring geodesics. 

Equation (A. 11) governs the stability of the geodesics and provides a convenient 
paradigm for describing the focusing and spreading of a geodesic bundle in M. The 
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tendency for neighboring geodesic to converge or diverge is determined by the matrix 
Kjj . If the curvature is positive definite in every plane containing the ray tangent, then 
neighboring rays from a common source point in all directions tend to converge and a 
three-dimensional ray bundle will focus. If the curvature is zero or negative everywhere 
along the ray, then neighboring rays from a common source point will tend to eventually 
diverge. (If Y 0 * 0 then the initial shape of the wavefront may cause the formation of a 
caustic independently of the focusing properties of the medium, as happens with light 
reflected from a concave surface in Euclidean space.) The conditions K n >0, and 
det Kjj > 0 along y(A) are necessary and sufficient for the quadratic form KjjYjYj to be 
positive definite. 

A.8 Conjugate Points 

Two points on a geodesic y(A) , P and Q, are said to be conjugate if there exists a 
nontrivial Jacobi vector along y(A) that vanishes at P and Q. According to the 
conjugate-point theorems of differential geometry, if upper and lower bounds exist such 
that k\ < R(Y,T,Y,T) < k 2 along some segment of the ray, then the period of affine 
parameter between consecutive conjugate points satisfies n/<T <n / y[k ^. If the 
sectional curvature is zero or negative everywhere along y F (A) , then neighboring rays 
will eventually diverge. (A single caustic may form if the initial region is concave but 
otherwise one does not expect the periodic development of focal points further 
downstream.) In situations where a ray occasionally may enter convergence and 
divergence zones, one cannot conclude from checking the Riemann tensor whether 
conjugate points will occur and how frequently. 
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Appendix B 

NULL GEODESICS AND ACOUSTIC RAYS 


This section introduces the connection between the null geodesics of a pseudo- 
Riemannian manifold and the bicharacteristics of the equations of hydrodynamics of 
ideal, isentropic fluids for the case of subsonic flow. The emergence of a geodesic 
structure from the study of hyperbolic partial differential equations is a hallmark of the 
method of characteristics. The same geodesic structure emerges in the study of acoustic 
rays, a connection that has been pointed out in the past by researchers in both the 
acoustics and relativity communities. In most cases one appeals to the high-frequency 
approximation although it should be pointed out that the same structure emerges in the 
general treatment of the system describing the propagation of discontinuities. The 
majority of this section contains little new information and is written for completeness, 
rigor and convenience, because the variety of different approaches has not been cataloged 
in one place. Every attempt at completeness is made while maintaining brevity. 

B.l The Bicharacteristics of the Equations of Hydrodynamics 

The equations of hydrodynamics for ideal, isentropic fluids consists of the continuity 
equation and Euler’s equation, respectively given by 


d t p + o-Vp + pc 2 V-o = 0 , 

(B.l a) 

d t v + u ■ Vu + — Vp + VO = 0, 

(B.lb) 


P 


where p is the fluid pressure, v the fluid velocity, p the fluid density, and where the 
fluid obeys an equation of state p = p{p ), with p' = c 2 defining the local sound speed. 
The potential <b is due to external forces such as gravity acting on the fluid. Equations 
(B.la) and (B.lb) are a system of first-order quasi-linear partial differential equations. 
Applying the method of characteristics to this system leads to the characteristic matrix 


D(p joc 2 {v<pf 

P~^<P 1 3x3 Dcp 


(B.2) 


where the function (p(x,t) is the wavefront, defined as a surface in space-time on which 
the initial data of the fields are specified. 
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The characteristic conditions are satisfied by setting the determinant of the matrix 
equal to zero leading to 

Q = {D<pf{{D<pf -c 2 V<p-V<p)=0. (B.3) 

Each term in Eq. (B.3) defines a partial Hamiltonian for the system and may be treated 
separately. Defining H l =(D<p ) 2 and H 2 =((j)<p) 2 -c 2 'V<p-V<p), one set of 
characteristics is determined by 

H 2 =g^d M <p8 v <p = 0, (B.4) 

where the contravariant metric tensor g“ v has been introduced with g 00 = -c 2 , 
g 0 ' = -ju t , g iJ = dy - Hi/J , , with ju = o/c defining the local Mach number, or in matrix 
form 


lf-1 

8 = 


c 2 8 n — v j Of 


The corresponding covariant metric tensor given by 


-c +o -v 

-a, 8,-f 


(B.5a) 


(B.5b) 


is the inverse of Eq. (B.5a) satisfying g fta g av = 8 V . A factor of pc 2 has been removed 
from Eq. (B.4) for simplicity but otherwise having no effect on the geometry of the 
curves defined by it. The bicharacteristics, which are identified with the curves along 
which energy is propagated, or with the rays of the system, are given by a set of 
trajectories in space-time. Defining an appropriate parameter X along these curves and 
introducing the canonical momentum p u = d^(p in H 2 , the bicharacteristic equations are 


dx M _ 1 dH 2 
dA 2 dp u 


(B.6a) 


dp„ _ 1 dH 2 

dX ~ 2 8x M 


\p.pp„g’ f - 


(B.6b) 


In Eq. (B.6b) the derivative of the contravariant metric can be replaced by the derivatives 
of the covariant metric by differentiating g fia g av = 8 l p , thus combining Eqs. (B.6a) and 
(B.6b) leads to the following equation for the bicharacteristics: 
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(B.7) 


d 2 x^ 
dX 2 


2 


g MV 


which is precisely the equation for the geodesics of a differentiable manifold M with 
metric g . From Eq. (B.7) and the characteristic condition, it follows that the 
bicharacteristic curves are equivalent to the null geodesics. 

B.2 Linear Acoustics 


To describe acoustics in the presence of background fluid motion, all fields in the 
problem are written in the form o = w + v, p = p 0 + p x , and p = p 0 + p x , where w , p 0 
and p 0 are respectively the velocity, pressure and density of the fluid medium, and v , 
p x and p x are contributions due to acoustic phenomena. In the following treatment of 
the problem it is assumed that the perturbations are weak, leading to a linear theory for 
the propagation of sound. No restriction is placed on the state equation for the 

background fluid motion, while the acoustic fluctuations are assumed to obey the 

adiabatic condition, p x = c~ 2 p x . The background fields are assumed to obey the 

equations of hydrodynamics; with body forces such as gravity lumped together in the 

background equations, independently of the acoustic field. Applying these assumptions 
to Eqs. (B.la) and (B.lb) gives 

{d*p x + c 2 p 0 V • v}+ {pjV • w + c 2 v • Vy o 0 + p x c 2 D*c~ 2 }=0 , (B.8a) 

{y o 0 D*v + Vp x }+ \c~ 2 p x D*w + y o 0 v ■ V>v}= 0, (B.8b) 

for the acoustic perturbations, where D* =d f +#-V, following Thompson’s notation 
and the adiabatic assumption has been used to eliminate derivatives of p x . It is not 
assumed that the background quantities or their variations are small. Letting the four- 
component vector u l>] denote the field variables with a [,] = v t for / = 2, 3, 4, and 
u m = p x , Eqs. (B.8a) and (B.8b) may be written in matrix form 


L (X) [u] + M u = 0, 


(B.9) 


where T (1) [---] is a first-order differential operator acting on u while M is an ordinary 
matrix acting on u. In both cases these matrices depend on the background fields and the 
acoustic fields: 


L 


(i) 


/) 

V 


hxiPoD*] 


M = 


c 2 V -w + D*c 2 
c 2 D*w 
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One can see by inspection that the first terms in Eqs. (B.8a) and (B.8b) resemble Eqs. 
(B.la) and (B.lb), leading to the conclusion that the characteristic condition for acoustic 
perturbations in the presence of background fluid motion is given by Eq. (B.5a), where 
the term s in the metric contain the background w instead of v . 

B.3 The Eikonal Approximation, Acoustic Rays 

Applying D to Eq. (B.la) after dividing through by pc 2 , taking the divergence of Eq. 
(B.lb), and taking the difference of the resulting equations, leads to 

D^-Dp-V—Vp + A 1= 0, (B.lOa) 

pc p 

while taking the gradient of Eq. (B.la) and applying D to Eq. (B.lb) leads to 

DpDu-Vpc 2 V-u + A 2 =0, (B.lOb) 

where A x = -d k v j d j v k and A 2 =-d k pVu k . Applying the same procedure to Eqs. 
(B.8a) and (B.8b) leads to a similar result for the linear acoustic field in the presence of 
background motion: 

\ D* — ^—D*p x - V • — Vpj i + = 0, (B.lla) 

l Po c Po J 

[d* p 0 D*v -Vp 0 c 2 V -v}+A 2 = 0, (B.llb) 

where 

Aj =D" UVV-W -V-+ Z>*(v-V In/?„)+£>’ \^-D*c~ 2 

l Po c J l Po c ) VPo ) 

- d k (v • S/w k )- 8 k VjdjW k , 

A 2 = D* ( p 0 v • V>v)+ D* {p x D*v>)~ v(pjV • wj- v(c 2 v • Vp 0 )- d k p x Vw k . 

The equation for the eikonal in the high-frequency limit of linear acoustics, like that for 
the bicharacteristics, comes from the highest-order derivatives. The terms Aj and A 2 
contain only the field variables and first-order derivatives. The eikonal approximation is 
defined by taking a—>0, p 1 =n a e" /,la and v = d , a e" /,la , with 7t a = ;r 0 + an x + • • •, 
a a =a 0 + aa l H—, where complex fields have been used for simplicity. Inserting these 
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expansions for the field variables and collecting all term s of order a 2 leads to the 
following eikonal equation in matrix form: 




(B.12) 


where H 2 is defined in Eq. (B.3) and B = c 2 Vcp <8>V<p - \ 3 x3 (d*< pf. Defining 
Q = -H 2 © B for the block-diagonal matrix in Eq. (B.12), the condition for the existence 
of a nontrivial solution to Eq. (B.12), det Q = -H 2 deti? = 0, leads to the condition 

Q = (D*<pY((D*<pJ-c 2 V<p-V<pf =0. (B.13) 

By inspection of the term s Aj and A 2 , one sees that they may be set to zero in the limit 
of a slowly varying environment since every term contains derivatives of the background 
fluid variables. This leads to a decoupling of the equations: 

D* — ^—D*p l - ■ —Vpj « 0, (B.14a) 

Po c Po 

D*p 0 D*v - Vp 0 c 2 V • v « 0. (B.14b) 


B.4 Time-Parameterized Rays in Three Dimensions 

The surface of constant time (p{x,t) t=t = (p(x ), generated by taking a constant time 
cross-section of the wavefront, is referred to as a phase surface (Fig. Bl). (From here on, 
a wavefront is a surface of equal time, and “in phase” is taken to mean equal time of 
flight.) 

In general, displacements in the null hypersurface are subject to the constraint 
g 00 dt 2 +g iJ dx‘dx i +2g 0i dtdx' = 0, which leads to the following constraint on dx/dt : 


gii dx dx J | ^ g 0 . dx 
goo dt dt g 00 dt 


(B.15) 


For the specific metric given in Eq. (B.5a), Eq. (B.15) states that the ray velocity relative 
to the moving fluid equals the local sound speed (dx/dt-o)-(dx/dt-o) = c 2 (see Fig. 


B2). 
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Figure Bl. Sample of a complete space-time ray structure. (A) Family of geodesics labeled 
by s, where the 5 parameterized flow lines are constructed from the deviation vector acting 
on a fiducial geodesic. (B) Each flow line represents a curve of constant affine parameter. A 
tangent vector and the deviation at a single point are illustrated. (C) A curve of equal-time 
(wavefront) constructed from the intersection of a constant time plane with the geodesic 
flow. Equal-parameter curves do not lie in the constant time surface in general. 



Figure B2. Relationship between the fluid velocity, 
ray tangent, and wavefront normal in Euclidean space. 
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An equation for dx/dt is derived from Eq. (B.7) by writing the equations for time and 
space separately: 


dx , dx‘ dx ] _ k dt dt „ k dx‘ dt . 

——+Tij -+r oo-+ 2r <0-= o, 

dA 2 dA dA dA dA dA dA 


d 2 t dx‘ dx j dt dt 0 dx' dt 

- t - + I ij - hi 00-1- 21 i 0 - = U . 

dA dA dA dA dA dA dA 


(B.16a) 

(B.16b) 


Changing the parameter along the ray from A to t leads to 


d 2 x k d 2 t/dA 2 dx k k dx‘ dx J k dx 1 

— —+ 7 —-— 7 T— + r k ij - +r oo + 2 r <o — = o, 

dt (dt/dXf dt dt dt dt 


(B.17a) 


d 2 t/dA 2 o dx 1 dx 1 , ~^o dx 1 

7-7T + i ij -1" t 00+21 io - = U . 

(i dt/dA ) dt dt dt 


(B.17b) 


Using Eq. (B.17b) to replace the second term in Eq. (B.17a) gives 


o ax ax ' o 0 ax \ax 

-1 ij -1-1 00+21 ,0 - - 


+ r* f ——+r*«.+2r *, 0 A = o. 

dt dt dt 


(B.18) 


Using the Christoffel symbols, Appendix D, and some algebra, the equation for the 3- 
dimensional ray paths as a function of time becomes 

x k = —— Sy (x ; - v t )(xy -Vj)+ (x k - v k jyd - lnc + 2f -Vine) 


+ d t v k -cd k c-<±> kj (x 7 [s kj -co kj \>j. 


(B.19) 


Equation (B.19) is converted to an equation for the wavefront normal through the 
identification cn = x-u , and can be used to determine either the evolution of x or n : 


—h k +Q 1H n i =0, 
dt k h 1 


(B.20) 


where Cl ki = ^£ jki (v + s xn + 2 Vc x n)j. 
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Along with a complete presentation of the bicharacteristics and a discussion of the variety 
of situations described by null geodesics, the main results of this appendix are Eqs. (B.7), 
(B19), and (B.20). In particular, note that the geodesics determined by Eq. (B.7) describe 
the propagation of discontinuities and the evolution of the wavefront in the most general 
situations, for nonlinear disturbances propagating in a random environment. The catch- 
22 is that in order to find these curves in the nonlinear case one needs the complete 
velocity field, u = v + w , which requires knowing the complete solution in the first place. 
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Appendix C 

CONSTRUCTION OF THE INTERNAL FRAME FIELD 


C.l Ray-Centered Frame in Space-Time 

To calculate the geometric spreading of an acoustic beam from the deviation equation 
the Riemann tensor must be projected into a ray-centered frame field. The frame being 
defined formally on the null hypersurface requires a basis of four mutually orthogonal 
vectors in space-time. Following the notation in Hawking et al? 2 , the ray tangent, T u , is 
chosen to serve as one of the coordinate directions. A second null vector, If , satisfying 
the condition IfT v g MV = -1, is chosen as the second basis vector. To complete the local 
geodesic coordinate basis, two space-like directions, e“ and ef, satisfying the conditions 
T a e Ia =0, L a e, a =0,and e“e Ja = 8 u for I,J= 1, 2, are introduced. 

The two null vectors define a time-like direction E t =c~ 1 ( 1 v) that is orthogonal to 
the two-dimensional space-like hypersurface (e ,) at all points in space-time. After 
defining this pseudo-orthonormal basis at an initial point on y F (A) , the basis is parallel 
transported along y F (A) to define a new basis at each point by solving the parallel 
transport equation: 

^A+rv7''(c,)' = o. (c.i) 

dA 

The null vector If is automatically parallel transported along the null geodesic . 

To ensure that neighboring rays are emitted from the source at the same “time” as the 
fiduciary ray, the initial space-like basis is chosen to be purely spatial in the coordinate 
frame, (e ,)" = (0 e/), and the initial deviation Y 0 M = (o Y 0 ) is chosen to lie along one 
of the initial basis vectors. 

Equation (C.l) formally preserves the constraints placed on the initial conditions. 
Beginning with 16 quantities, 4 vectors with 4 components each, the conditions impose 
11 constraints leaving 5 independent quantities. Four of these are the components of T M 
determined by the geodesic equation, leaving only one remaining quantity to calculate. 
The vector If may be determined algebraically in terms of the environmental 
parameters. The remaining degree of freedom may be identified with a two-dimensional 
rotation of the space-like basis {e ,} in the (e ,) hyperplane expressed as Rjje J0 . 
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C.2 Ray-Centered Frame in Cartesian Space, Wavefront Normal Basis 


In this section an alternate choice of internal basis is presented. This representation, 
termed auxiliary basis, is better suited for visualization and interpretation of the 
deviation, as it always remains tangent to the surface of constant phase in space. From 
the conditions T M e v I g^ = 0 and e^ejg fiv = 1, one may verify that 

n-(e,-e° I t)= 0, jjitf - te° || = 1, (C.2) 

with respect to the ordinary dot product in three-dimensional space. Hence the set of 
three-dimensional vectors 

~o dx k 

e i = e i ~ e i — (C.3) 

at 

is tangent to the wavefront everywhere along the ray and the basis (n, e r ) remain s 
orthonormal as it is propagated along the ray. The vectors e, are referred to here as an 
auxiliary basis. An equal-time deviation is constructed by projecting the components Y t 
onto the auxiliary basis. The deviation projected into the different basis will point in 
different directions in four dimensions but the magnitude is independent of the choice of 
basis, Jr| = ^ Y m Y m = ^/ij 2 +Y 2 2 . 

The sectional curvature needed for the spreading equation may be rewritten as 
4 K u = R flvaf A, M 'd \, where A/' 1 ' = T M e v r - T v . The mixed components of this 
tensor, 

A j 0k = T°e k -T k e°j = T°(e k -e] (< ix k /dt )) = T°e k , (C.4a) 

are clearly proportional to the auxiliary basis. The pure space components, 

A/ k = T'e) -T k e k = Te k - T k e ], (C.4b) 

are equal to the components of the cross product f x in three dimensions. The tensor 
A/ v is parallel transported along y{X ), 

—A/ r + rV" A/ v + Y v a pT a A I Mp = 0 . (C.5) 

dA, 

Equation (C.5) tracks six fields, counting indices in four dimensions. It is clear from the 
definition of A/' 1 ’ that there are only three degrees of freedom. By separating Eq. (C.5) 
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into individual components, the following transport equation for the auxiliary basis e/ in 
E 3 is derived: 


— e? +-£ ih (vxL ) + sxn + 2Vcx n) el = 0, 
dt 1 2 jk,K hi’ 


(C.6) 


with s k = HjSjb . Note that an overall factor of dt/dA has been divided out. In principle, 
the parallel transport equation has a solution in terms of a path-ordered integral involving 
the SO(3) generator evaluated along the ray. In practice, these are not known in closed 
form. However, assuming that the ray equation has been solved either analytically or 
numerically the normal n is known, reducing the problem to a two-dimensional rotation 
about the wavefront normal followed by an application two Euler matrices to orient the 
basis along the ray path. 

Choosing the initial basis to match a global Cartesian-coordinate basis 
(n(0) e ] (0) e 2 (0)) = (i j k) , the rotation matrix taking n( 0) —» n may be found 
explicitly. The rotation of (e)) about n is denoted R(a) . The action of this sequence 
on the initial basis vectors, interpreted as an active transformation on the initial values 
(n(0) e, (0)), is denoted U = R((p)R{6)R(a) . The rotation R{a ) has a block-diagonal 
form R(a) = 1 ® r , where r is a two-dimensional rotation about the x-axis of the global 
coordinate system. Inserting this into Eq. (C.6) leads to the following differential 
equation for r: 


(C.7) 


where 


■C 


The solution to Eq. (C.7) is 

fcos a -sinci^ 
l sin a cos a J 


(C.8) 


where a = j co\^ dt , with the integrand evaluated along the specific ray path in space. The 

dependence of a> on specific Cartesian indices comes from the choice of lining up the 
initial basis with the global coordinate system. 

Equation (C.8) reduces the problem of tracking eight degrees of freedom plus 
constraints to one of evaluating a single integral. This reduction is not without cost. It is 
clear that the denominator will vanish whenever a ray has a horizontal turning point. 
While this may not happen often in modeling underwater systems where attention is 
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focused on parabolic rays trapped in sound ducts, it will happen in more general 
situations. The danger of this divergence occurring indicates that it may be more 
practical to track the components of the larger dynamical system, using the constraints to 
regulate error along the way. 

The ability to reduce Eq. (C.l) to Eq. (C.7) depends on the identification of A ^ with 
the wavefront normal basis. In general, this trick may always be employed to reduce the 
number of quantities needed to calculate the sectional curvature and is invaluable in 2- 
dimensional ray systems, 3-dimensional space-time systems, as it allows one to generate 
e by rotating h by 90 . Without this insight one would be compelled to either track all 
three components of e by Eq. (C.l) or evaluate Eq. (C.8). 
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Appendix D 


RIEMANN TENSOR AND CHRISTOFFEL SYMBOLS FOR THE ACOUSTIC 
METRIC 


For the acoustic metric defined by the characteristic condition, Eq. (B.5), 


g MV 

1 r-1 -V t \ f-c 2 +o 2 -u '| 

c 2 (j-Uj c 2 S ji -v j v i y Sftv y — Vj S Jt y 


direct substitution into Eqs. (A.4) and (A. 10) lead to the following: 

r° 00 = d 0 \nc + ^-S, , 

2c 2 9 

(D.la) 

r °, 0 

= d:\nc--~fS* & 

' 2c 2 9 

(D.lb) 

r ! 00 

= vd Inc - dp, + cdp-—u m, - J -r (c 2 8 kj - v k v, )s ik , 

1 ~ v ‘ 1 1 2 1 11 2c lX kl k l> Jk 

(D.lc) 

r'yo 

_ , 1 vp k 0 

= up. In c + — op — '-Pr S ik , 

11 2 9 2c 2 Jk 

(D.ld) 

r % 

and 

= u k r° i j=-^ T s ii , 

2c 

(D.le) 

R„ imj 

=^(s..s, ) -s M s,) , 

(D.2a) 

R mj 

= toc-Std, \nc)+-^(s lt s J ,-s l s h ,) 

, (D.2b) 


R 0i0i = ~-d t S, ~-v n (d i S in + 5 ,op, )+ cd t 8 ,c + -S ii d t In c -°^S w S im 
u '"J 2 ' 1 2 v 1 1 1 2 4c 2 1 

+ L j( S jn di Xac + S in d ] lnc - s ,j d n +(o im S jm + co Jm S im ) . 

(D.2c) 
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Typically the Riemann tensor contains second derivatives of the metric components as 
well as first derivatives squared. A significant feature of this system is the lack of 
second-order time derivatives of the environmental parameters and the coupling of 8 t c to 
the fluid velocity. For a stationary medium the focusing is determined solely by the 
spatial derivatives of the sound speed. 
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Appendix E 

FERMAT’S PRINCIPLE AND SNELL’S LAW 


E.l Fermat’s Principle 

Fermat’s principle follows immediately from the line element by solving for dt : 
yj(v-df ) 2 +(c 2 -v 2 )df -dr -O- dr 


(E.l) 


where the + root is taken for subsonic flow. For supersonic flow, dr ■ v > 0 and the 
travel times are determined by the - root. 

E.2 Metric Symmetry and Snell’s Law 

An isometry (symmetry of the metric) is indicated by a cyclic coordinate 30 . The 
presence of a cyclic coordinate in the metric, labeled x c , leads to a conservation law for 
the corresponding momentum or conjugate variable, 


dp c dj^ £ r>> 


dX 


d )l 


dX 


(E.2) 


The vector called a Killing vector, points in the direction of the x c coordinate 
curves. (The statement “...presence of a cyclic coordinate” x c is equivalent to the 
absence of that coordinate from the metric altogether.) When there is enough symmetry, 
the second-order geodesic equations are reducible to a set of first-order equations. 

Specializing to a time-independent environment, the Killing vector = (l 0) leads 
to the conservation law 


( 2 2 \<*t _ dr 

-l c -v ) - v -= /f n . 

V ’dX dX 0 


(E.3) 


Traditionally, this is associated with the energy of the particle, and for the metric 
signature in Eq. (B.5) /c 0 is negative. (From here on, /r 0 > 0 and an overall minus sign is 
dropped from Eq. (E.3).) Consider the problem of a layered medium where c and o 
depend on only one Cartesian coordinate, z. There are two more conserved currents, 
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which may be thought of as components of the ray’s translational momentum. The 
Killing vectors i and j lead to 


dr dt _ 

-1> 7 

dX dX 3 


= K, 


(E.4) 


where k = K t i + ic 2 j and o T = o x i + o j . Imposing the null constraint on the components 
of the tangent vector leads to 

-c 2 i 2 +K-K + (z-o z tf = 0 (E.5) 

and 

cH = ic 0 -v t -k—v z (z-vJ). (E.6) 

Equations (E.3) and (E.4) may be used to derive a generalized version of Snell’s law 
originally presented by Kornhauser 14 for the special case v T = vi . By definition Eq. 
(E.3) states that k = tch T which in turn states that the projection of the wavefront normal 
in the x-y plane points in a fixed direction. Equation (E.3) may be rewritten to give 
k 0 =f-ch = tc(c + v ■ h) , relating the energy to the angle between the wavefront normal 
and the spatial component of the ray. When o =o T = ui, Eq. (E.6) becomes 
k 0 = ic(c + on x ). Taking the momentum-to-energy ratio in this case gives Snell’s law for 
the wavefront normal: 


cos# 

c + t>cos# 


(E.7) 


with a = K x /K 0 and n x =cos6. The ray angle is related to 6 by ccosd + u = rcosy/ . 
The magnitude of the ray velocity follows from the null constraint 

r = ucosi// + -yjc 2 -l> 2 cos 2 ^ . (E.8) 


Insertion of Eq. (E.8) into Eq. (E.7) gives Kornhauser’s result for Snell’s law: 

_1 -s sin 2 y/+ cosy/yll -s 2 sin 2 y/ 
c 1 - £ 2 sin 2 y/ + cos^-y/l-^ 2 sin 2 y/ 


(E.9) 


The paths of the rays propagating in a layered media with o, = 0 are 
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(E.10) 


dt _ (/r 0 -u t -k) 
dA c 2 

^- = k + ^(k 0 -u t -k), (E.ll) 

dA c 


dz_ = ± ^(k 0 -v t -k) 2 -k-kc 2 
dA c 

Taking appropriate ratios of Eqs. (E.10) and (Ell) leads immediately to the ray-path 
integrals encountered in underwater acoustics. 

The procedure may be generalized to any set of coordinates. A derivation of two- 
dimensional Snell’s law in polar coordinates is now presented. The acoustic line element 
in polar coordinates is 

-(c 2 -o 2 )dt 2 - 2o r drdt - 2v e rd6dt + dr 2 + r 2 d0 2 = 0. (E.13) 

Specializing to a rotationally symmetric, time-independent environment with v =v0, 
leads to the following set of equations for the ray: 


c 2 r 2 r 2 = -a 2 + (r - av ) 2 , 

(E.14) 

c 2 r 2 0 = ac 2 + v(r - av ), 

(E.15) 

c 2 rt = (r- av) . 

(E.16) 


Following all the same steps leads to the same result for Snell’s law, Eq. (E.9), with 0 
interpreted as the instantaneous angle between the wavefront normal and 0. 
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